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Abstract 

We describe a general Godunov-type splitting for numerical simu¬ 
lations of the Fisher-Kolmogorov-Petrovski-Piskunov growth and dif¬ 
fusion equation in two spatial dimensions. In particular, the method 
is appropriate for modeling population growth and dispersal on a ter¬ 
restrial map. The procedure is semi-implicit, hence quite stable, and 
approximately 0(Ax 2 ) + 0(At 2 ) accurate, excluding boundary con¬ 
dition complications. It also has low memory requirements and shows 
good performance. We illustrate an application of this solver: global 
human dispersal in the late Pleistocene, modeled via growth and dif¬ 
fusion over geographical maps of paleovegetation and paleoclimate. 

Keywords: reaction-diffusion equations , finite-difference solvers, 
population dynamics 


1 Introduction 


There is an increasing interest in modeling population dynamics at large 
spatial and temporal scales, for example the modern human out-of-Africa 


dispersal ( 

Eriksson et ah, 

2012; 

Henn et al., 

2012; 

Nikitas and Nikita, 

2005; 

Young and Bettinger, 

199, 

5) or 

Neanderthal dispersal and extinction ( 

Calle- 


gari et ah, 2013). These models are required to interpret local and global 


patterns of genetic, phenetic and cultural variation (Bouckaert et ah, 2012 
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2004). 


Fisher (Fisher, 1937) studied the description - via a reaction-diffusion 


equation - of an analogous but one-dimensional problem: the propagation of 
an advantageous genetic mutation within an already-present population, sit¬ 


uated along a coast line. Kolmogorov, Petrovskii and Piskunov (Kolmogorov 


et ah, 1937) were more general; in particular, their analysis treated the two¬ 


dimensional case. Such a model (called Fisher/KPP in the following) was 
first applied to the dispersal and growth of a population by Skellam (Skel- 


lam, 1951), and serves as an important control for designing and validat¬ 


ing other more complex spatiotemporal population models (Callegari et al. 


2013). Coupling population dynamics with models of large-scale changes 


in continental topography, climate, and ecosystem productivity is essential 
to understand the role of environmental constraints on patterns of genetic, 


phenetic, and cultural variation among human populations (Callegari et al. 


2013). 


Here we present a stable and efficient finite-difference solver for the Fisher/ 
KPP equation on 2-D domains of arbitrary shape (e.g. geographical maps), 
and show how it can be extended to include environmental fluctuations. In a 
brief outline of our paper, we will: review the derivation of the Fisher/KPP 
equation (Section [2]); develop finite-difference schemes in 1 and 2 dimen¬ 
sions for constant environmental carrying capacity, K, (Section |3j); extend 
the scheme to allow for space- and time-dependent /C(x, t) (Section |4j) and 
irregular domains such as geographical maps; show an application of this 
technique to the out-of-Africa dispersal of Homo sapiens by using net primal 
productivity (NPP) as a proxy for /C(x, t) (Section [5]). 


2 The Fisher/KPP equation 


An intuitive way to get the Fisher/KPP equation (Fisher, 1937; |Kolmogorov 


|et al.[[l937 ; Skellam, 1951| ) is as follows. A current j of particles (e.g., individ¬ 
uals) moving across an interface located at x is proportional to the gradient 
of the population density p (“Fickian diffusion”) 


j = —c Vp. 
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The rate of change of p is then given by the mass balance equation (Reichl 
2004), which for Fickian diffusion reads 


dp 

dt 


— c V 2 p = p. 


If p = 0, this is the heat equation when c = 22/2 and V is the diffusion 
coefficient. For lack of a better model, we assume c is a constant (Young & 


j^)p + cV 2 p, 


( 1 ) 


Bettinger (Young and Bettinger, 1995)). The source term p is usually mod¬ 
eled by a logistic growth function, p = \p(l — p/K ), yielding the Fisher/KPP 
equation 

dp _ . (, p 
dt V 

where K, is called carrying capacity and A is the growth rate. In (|TJ) we 
assumed JC is constant, but it suffices that there is an upper limit /C, in 
which case 0 < p/iC < 1. This scaled version will be used in ([ 2 ]) below. Fisher 
and KPP were particularly interested in the traveling wave case, p(t , x) = 
/(x 0 + vt). Notice what happens here if / exists: 


v • V/ = A1 


{)/ + cV 2 /, 


which in 1-D becomes a second order ordinary differential equation 

if . (, A, dV 

”s = H 

If v were known, this ODE could be solved using pvp4c from MatLab, for 
example. By the rescalings show in Table [lj for constant K, the Fisher/KPP 
equation Q will be used in the form 


du . 1 9 

¥ = (l-«)« + -Vu, 


( 2 ) 


where the only sensible solutions have 0 < u < 1. The initial distribution 
w(0,x) = Uo( x ) must be defined for all x. 

In Murray (([Murray, 2002), eq. (11.17)) our growth coefficient A is called r 
and c is denoted by D, whereas in Young and Bettinger (|Young and Bettinger 


1995) the growth coefficient is R and the diffusion coefficient is K. These 
inputs to our code are given in units of yr -1 and km 2 /yr respectively. 
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Table 1: Variables in left column are scaled versions of those in right column. 


In 

2) 

In (] 

L) 

X 

t 

u 

vT x 

A t 

p/JC 


3 Numerical methods and splitting 

In one dimension, ([2]) can be solved using the MatLab function pdepe. In 
fact, if the system is two-dimensional but rotationally symmetric, pdepe can 
again be used with the radial part of the Laplace operator in cylindrical 
coordinates, 

o 19 0 . 

V = r——b non-contributmg terms, 

r Or or 

requiring only that one sets a pdepe parameter m=l. Although the MatLab 
function pdepe is robust, it cannot be generalized to arbitrary 2-D domains. 
However, it is a valuable control for testing more general solvers, and a more 
general solver is what we wish to explore here. 


3.1 The finite-difference scheme 


Since the map on which we will be working is a pixelized plane, an obvious 
method uses finite differences. First, however, let us examine the 1-D case 
for ([2]). In this situation, the second order derivative becomes a differencing 
operator in matrix form acting on the vector {uj, j = 1, n}, where Uj = 
u{x 0 + (j - l)Ax), 

d 2 u 1 

_ _V _ An, 

dx 2 (Ax) 2 ’ 

where the matrix A is 


A = 


-2 
1 
0 
0 
0 


1 0 

-2 1 

1 -2 

0 1 

0 ... 


0 

0 

1 

0 


1 


0 \ 
0 


- 2 / 


(3) 
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If h is the time step, the Courant-Friedrichs-Lewy (CFL) parameter (Strang 


and Aarikka, 1986) is 


k = 


h 


2 (Ax) 


An explicit integrator for ([ 2 ]) would require k <1/4 (LeVeque, 2007; Strang 


and Aarikka, 1986). In our case, because the boundary conditions are so 


irregular on a map, we are less interested in a method of higher order than 
2nd because stability is more important (Godunov and Ryabenkii, 1987). 

Using this notation, the lowest order approximation is Euler’s method 
which estimates the next step u(t + h) by 


ue = u(t) + k Au(t ) + h (1 — u(t))u(t), (4) 

which should be considered a vector equation in u(t) = {uj(t),j = 1, n}. The 
logistic terms, which are diagonal, should be taken to mean ((1 — u)u)j = 
(1 — Uj)uj for j = 1 . 2 ...., n. Euler’s method is both low-accuracy and 
usually unstable if it is used alone over many steps. But, it is 0(h) accurate 
and thus useful as an explicit estimate in 0(h) terms. An application of the 
trapezoidal rule yields 

k 

u(t + h) = u(t) + — (Au(t T h) T Au(t)) (5) 

Zj 

+~ ((1 - u(t + h))u(t + h) + (1 - u(t))u(t)) 

and is an 0(h 2 ) + 0((Ax) 2 ) accurate procedure but solving the quadratic 
vector equation ([ 5 ]) for u(t + h) is awkward. To the same 0(h 2 ), we propose 
a semi-implicit procedure which uses the Euler estimate 0 to replace one of 
the terms in ([5]): 

k 

u(t + h) — u(t) H— (Au(t + h) + Au(t)) (6) 

2 

h 

+- ((1 - u E )u(t + h) + (1 - u(t))u(t)). 

Equation (|6]) can now be solved as a linear system, 

(l - \ A ~ ~ u(t + h) = u(t) + ^ Au(t) + ^(1 - u(t))u(t), (7) 
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because the matrix, 1 — |A — |(1 — ue), on the left hand side is explicit, 
as is the right hand side. That is, this matrix and the right hand side 
contain only old data, namely only information from the previous step, u(t). 
Euler estimate ue is an explicit one step computation using u(t). Significant 
advantages are: the matrix on the left hand side is tridiagonal with constants 
on the sub/super-diagonals, and the diagonal terms are 0(1) strong. The 
procedure 0 is only linearly stable but we will show empirically that it 
gives good results when compared to pdepe when this MatLab function is 
appropriate, that is, in both the one-dimensional and rotationally symmetric 
2-D case. Not only is the method 0 step-wise stable but also stable for 
initial data which may not be smooth. 

Figure [I] shows the results for h = 1/5, k = 2.5 compared to pdepe. 
Notice that at t = 20 the agreement is remarkable; and that at t = h, 
where the wave front profile is very steep, our Godunov splitting described 
in Section 3.2, specifically eq. (|8j) , is very stable. The CFL number, k = 2.5, 
used to get Figure [l] is much larger than would be possible with an explicit 
method (LeVeque, 2007). 


3.2 2-D case: Godunov—Strang—Yoshida splittings 

It turns out that a generalization to the 2-D problem is a straightforward 


variant of Strang-Yoshida splittings ( 

Strang 

1968; 

Yoshida, 1990), which are 

themselves variants of Godunov’s method ( 

Godunov and Ryabenkii, 

1987 

)• 


The following is a fully implicit variant of our two-dimensional scheme, with 
two intermediate arrays, u* and u**, 


k 

u* = u(t) + — (. A x u * + A x u(t )) (8a) 

u** = u* + ^ ( A y u ** + A y u *) (8b) 

+ ^ ((i - u**)u** + (l - «*K) 

k 

u(t + h) = u** + - ( A x u(t + h) + A x u**). (8c) 

In ([8]), the operators A x and A y are the same as ^ for directions x and 
y, respectively. For simulations on a lattice, up (t) = u(t, xq + (i — l)Ax, yo + 
(j — 1)A y), where 1 < i < N x , 1 < j < N y and Ax = Ay, the following gives 
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the action of the A x , A y operators: 


A x u hJ Ui—ij '2ii l : j T u i+ \, 

■Ay'U'iJ I Tpj “I - . 

The compression scheme and code outline given in Appendix [A] show that 
only a maximum of one row or column (i.e., max(A^, N y )) of storage is needed 
for u * and u**. 


Again because the fully implicit quadratic vector equation in (8b) is awk¬ 


ward to solve, we use an Euler estimate in one of the terms. Here is one 
integration time step of (j8j) in discrete semi-implicit form: 


1 - u* = ^1 + u(t) 

ue = u* + kA y u* + h{ 1 — u*)u* 


1 - ~Ay - -(1 - «** - fi + ~Ay + -(1 - U 


^ ' u * 


1 — — u(t + h) — ^1 + — A, 


u 


(9a) 

(9b) 

(9c) 

(9d) 


Equations (9a), (9c), and (9d) are solved in sequence as multiple independent 


tridiagonal systems for «*, u** and the final step u(t + h). 


3.3 Symmetries in 2-D case 

Our Godunov scheme i§ is not rotationally symmetric, and thus one way to 
estimate the error is to assess a solution using ([9j) for a symmetric problem. 
Again, we can use pdepe but now with the cylindrically symmetric parame¬ 
ter choice m=l (see Section 12.5 in (Higham and Higham, 2005)). Figure [2] 
shows that any asymmetries are not apparent without more careful exami¬ 
nation. Even the wave front portrait of the 2-D case in the left-hand panel 
of Figure [3] and the error estimate in the right-hand panel of the same Figure 
are not sufficiently quantitative. In particular, there should be no distinction 
between x and y directions in (|8]), while a 2-D plot of the error distribution 
shows a small asymmetry (compare the right-hand plot in Figure [4] to the 
left). 

For this reason, we implemented an alternating direction method, a la 


Crank-Nicholson (LeVeque, 2007; Ritchmyer and Morton, 1967), which makes 
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the error distribution more symmetric. The left panel of Figure [4] shows 
that, while symmetrization only slightly improves the r.m.s. and maximum 
errors, they have now the desirable property of being more rotationally equi- 
distributed: respectively, cyclic groups C '4 vs. C%. Relatively larger devia¬ 
tions from the pdepe solution now correspond to directions diagonal to the 
spatial lattice, as expected, and do not reflect the arbitrary choice of x and 
y in the integration. 


4 Fisher/KPP on maps 

The next natural step when applying a reaction-diffusion equation to the 
modeling of population dispersal is to include geographical and environmen¬ 
tal effects. In this Section, we discuss how to implement our solver on do¬ 
mains with space- and time-dependent 1C, and then how to treat irregular 
boundaries that arise when solving Fisher/KPP on geographical maps. Using 
the same Godunov-type splitting described above, it is more straightforward 
to do the simulations on a map than might be expected. 


4.1 Maps with space-dependent capacity 

In our scaling of ([2]), the maximum population density at x is unity. Thus, 
in the following we will use a scaled carrying capacity 0 < /C(x) < 1. Our 
algorithm ([8j) can be modified in a straightforward way for the case that 
/C(x, t) also depends explicitly on time: see Section 4.4 First, however, let 
us deal with the time independent /C(x) case, 


du 

dt 


(1 


u x 

K )U + 



( 10 ) 
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(we dropped the x dependence of K for convenience of notation), for which 
the Godunov splitting (|8j) now becomes 


1 — — Ar'j u* — ^1 + — A x ) u(t) 


u 


u E — u + kAyii + h(\ — ) , u 

/C 

, k , u E ,\ ,, /, k . h. u*.\ , 

= ( 1+ ^h (i -cr 

1 — — u{t + h) = (l + —ix**. 


(lla) 

(lib) 

(11c) 

(lid) 


Again, as in ([8]), multiple tridiagonal system solves must be carried out: N y x- 
direction solutions (lla), N x ^-direction solutions (11c), and finally, another 
N y ^-direction solutions (lid). Hopefully no confusion will result from the 
notation: k is the CFL parameter, while 1C is the (space-dependent) carrying 
capacity. 


4.2 A desert test of space-dependent capacity 

Now we are in uncharted territory. To assess if the solver © works, let us 
examine a problem where we can compute a solution by independent means. 
The test setup follows below. Its motivations will be explained further in 
Section 14.31 
The desert test: 

1. for — y 0 < y < y 0 , K- = K,(x) is independent of y\ 

2. for xq < x < xl-, let 1C — 1, while for xy < x < xh, let 1C = f r , where 
the fraction 0 < f r < 1 defines a desert (inhospitable region) in the 
domain. Finally, for xh < x < x i, again set 1C = 1; 

3. initialize u(x, y.t = 0) to a strip midway between x 0 and xl, then run 
the simulation to study the traveling wave behavior across the [xl, xh] 
desert. 

In other words, the variable carrying capacity domain has K = 1 for 
x < xl and x > xh, but 1C = f r in an x-direction desert. If the initial data 
w(x, 0) are widely distributed enough (nearly full y- width), near the middle 
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of the domain, i.e. y = 0, the problem looks basically one-dimensional. Thus 
we can again use pdepe from MatLab to compute the behavior of the one- 
directional wave as it passes through the desert, and compare this solution to 
the behavior of our Godunov method near this same center line. According 


to KPP (Kolmogorov et al., 1937), the velocity in the desert is the same as in 


the /C = 1 region, and far enough from transients it should be approximately 

V = V2a, 

where parameter a = ^(it = 0). In our case p{u) = u(l — u/IC), so ^(0) = 1. 
Thus, when not entering or leaving the desert, the velocity should be a/2 
(Kolmogorov et ah, [1937). Figure [5] shows that the velocity agrees with the 
KPP prediction, and is nearly constant except for short transients entering 
and exiting the desert, as expected. Our metric for measuring this velocity 
is to find, on the leading edge of the wave front, the position aq /2 where 
u(x i /2 ) = /C/2. 

If the jump in K, is deep enough (f r <C 1), however, the integrator will 
fail without some regularization. The parameters for the results shown in 
Figure [5] are /C = 1 above and below the desert, but /C = 0.01 in the desert. 
For this case, / | gives an instability with unpleasant sign changes, and 

a regularization scheme has to be used. We show this in the next subsection. 


4.3 Regularization against holes and noise 

An obvious problem with deep holes, jumps or ragged noise in the carrying 
capacity K. is this: the right hand side of equation (TlcJ) (as well as (16c) in 
the next Section) contains the term 


k . h._, 

1 +2 >1 » + 2 (1 -iC ) 


u 


( 12 ) 


which for large u *//C basically determines the sign of the u(t + h ) on the left 


side of (11c) (likewise (16c)). Since both u* and /C are positive, if /C is very 
small in some pixel, then we may have 


u* 2 

K > V 


in which case u(t + h) changes sign. This is unphysical, so we would like to 
regularize the term (|12[) . To do so, we need to find a monotonically increasing 
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function, call it g(x), such that 


g(x) ~ 


x when x is small 
1 when x is large. 


Multiple choices are available, as shown in Figure [6} We want a straight line 
with slope 1 when x is small, then a smooth but rapid cut-off when x = ^ 
gets too large. Some suitable choices are (1 — exp (—x 13 )) 1 ^, or just as cheap 
to compute, 


g(x) = (tanh (x^)) 1 ^ , 


(13) 


or 


any variant shown in Figure 6| We choose g(x) with j3 — 4. The regu¬ 
larization (13) to be used in (12) and thus ( |llc[ ) (likewise ( |16c ) in the next 
section) is then modified to 


£ - >£> (l4) 

where h is the step-size. In the desert test presented above, using the regu- 
larizer yields the same results as decreasing h tenfold. 

One alternative to the above regularization 0 is to use smoothing by 
a low pass filter which weights a center pixel (map coordinate i,j ) and its 
nearest (2 L — 1) x (2 L — 1) — 1 neighbors. A neighboring pixel with X, Y 
distances i x , j y from has weight 

HixJy) = (! - (ix/L) 2 ) (l - (. j y /L ) 2 ) (15) 


for all — L < i x < L and — L < j y < L , including center at i x = j y = 0. 
Neighboring pixels with coordinates (i + i x , j + j x ) having zero value (e.g. 
water), i + i x > N x or * + i x < 1, are ignored. In Mercator projection maps, 
Y coordinates are periodic in the longitude direction). For each (i. j)—pixel 
to be smoothed, a total of each accepted (non-zero) neighbor’s weight was 
kept and the resulting total was normalized appropriately. The choice (15) 
is an approximate Gaussian weight exp (—Ax 2 — Ay 2 ), cut off at distances 
| i x | > L and \j y \ > L or at uninhabitable pixels. Figure |7| sh ows an example 
of the effect of a low-pass filter (5 cells half-width, eq. (|l5|)) smoothing on 
one of the maps used in this study. A comparison between regularization 
and smoothing is shown on the right hand plot of Figure [9} 
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4.4 Time-dependent capacity maps 

Fluctuations in climate produce environmental changes in vegatation, sea 
levels, opening/closing of land bridges, waxing/waning of ice sheets, and per¬ 
turbations to habitable areas in general. Thus, time-dependent environments 
compel us to extend our procedure 0 for both space- and time-dependent 
/C(x, t) (see Section 5.2). 

Since 0 is basically the trapezoidal method (see section 5.3 in (LeVeque 
2007])), the modification for a time-dependent 1C is as follows: 


1 “ u * ~ f 1 + 4 A - 


u(t) 

ue = u* + kA y u * + h( 1 


u 


m 


)«* 


(16a) 

(16b) 


1 _ ] 1 A h d u e 
2 y 2 JC(t + h) 


))«**— ( 1 + -A y + -(1 - ^-y) ) u * (16°) 


1 — — u (t + h) — ) u **, 


(16d) 


where we have again suppressed the x dependence of /C(x, t) for simplicity 
of notation. 


4.5 Fisher/KPP on geographical maps 

In order to solve Fisher/KPP on an irregular domain such as a geographical 
map, it is sufficient to break down the map into x- and ^-direction segments, 
imposing aa = 0 boundary condition at their endpoints; the solver can act 
on each segment independently, alternating the direction (jLeVeque, |2007t 
Ritchmyer and Morton, 1967) of integration as discussed above. This ap¬ 


proach lends itself also to efficient parallelization. Appendix [A] and Figure [12] 
illustrate in detail our scheme with a sample Matlab code. 

Note that holes in K maps can represent the partial closing of land bridges 
without the necessity of re-segmenting land portions, as in Section [4] and 
shown in Figure [12] However, the regularization scheme (Section |4.3| ) should 
be used with caution. Att = 0 boundary condition region is not the same 
as one with low 1C, which can pass a tiny population into a subequent region 
where it may well flourish. For example, historically known falling sea levels 
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opened passages across Bab al Mandeb to open South Asia for human dis¬ 
persal (Siddall et al., 2003), and the retreat of the North American ice sheet 
opened a passage on the Bering strait around 12 kya (kilo-years ago). 


5 World-wide hominin dispersal 

We now turn to a sample application of the methods presented above: the 
world-wide dispersal of Homo sapiens out of Africa. 


5.1 Capacity maps 

Our construction of a time-dependent 1C uses Net Primary Productivity 


(NPP) as a proxy (Eriksson et al., 2012). The Miami model (Grieser et al. 


2006) was originally formulated in 1972 to estimate NPP (in grams (of car¬ 
bon) in dry organic matter/m 2 /day) from annual temperature and rainfall 
(Lieth et al., 1975). In order to compute our NPP maps, we obtained the 
temperature and precipitation data from simulations by the bridge program 


(bridgeproj 

) organized at the University of Bristol ( 

Bigelow et al.[ 

2003 


Harrison et al., 

2001; 

Pickett et al., 

2004; 

Prentice and Jolly, 

to 

o 

o 

o 


The 

simulation data that we used were computed on a 96 x 73 grid, which we 
interpolated to size 100 x 50 and converted to NPP maps by applying the 
formulas given in (Grieser et al., 2006). In Figure [7j the original NPP units 
in (grams of C)/m 2 /day were adopted. 

World-wide NPP data are difficult to obtain, so our Miami model-like 


maps are rough. As we showed in Section |4.2[ our Godunov solver is fairly 
robust with respect to abrupt steps in the carrying capacity /C(x). However, 
an alternative to the regularization scheme proposed in Section |4.3| would 
be desirable, when dealing with maps in which noise and holes may not 
necessarily correspond to real physical conditions of the model system. 


5.2 Time interpolation of maps 

We assembled 61 NPP maps, from 120 kya to 1 kya. These are in 4 ky steps 
for the first 10, 2 ky steps for the next 21, then 1 ky for the remainder. Since 
the time stepper in our Godunov scheme can have no information about the 
continuity in time of NPP maps, an interpolation scheme needs to be used. 
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If some estimates or proxies for the carrying capacities ICl, JCh at times 
tij tn are available, one possible estimate for )C(t) is a homotopy 

JC(x,t) = /C L (x)( 1 - S(t)) + JC H (x)S(t) 

where a sigmoid function 0 = S(ti,) < S(t ) < S(tn) = 1 will smoothly 
interpolate between the two time frames. There are many choices available, 
such as that used in (Eriksson et ah, 2012). For this study, we used the 
following variant. 

Start with the classical sigmoid 


S(j) = (l+e-‘) 


-1 


(17) 


which is zero at z = —oo and unity at z = +oo. The —oo < z < +oo 
interval is not what we want, but with a small change the following t z 
transformation permits many variants: 


z = 


2A T(t - t L ) - AT 2 
((* - (t - t L ))y 


(18) 


where AT = £# — Notice that S(z(ti, )) = 0 and S(z(tn )) = 1. The 


exponent u in (18) gives some freedom in choosing a particular form for z 
for almost any u > 0. If u < 1/2, d 2 S'/dt 2 has more than two sign changes, 
so v > 1/2 is preferable. By the choice v = 1/2, the interpolant is nearly 
a straight line: see Figure |8j However, its turn-up at t = tL and turn-down 
at t = tn numerically resemble very quick derivative changes. Thus, for this 
test case we choose u = 1. At both ends, all derivatives of S(z(t )) in t vanish. 
Also notice the forward/backward symmetry S(z(tn — £)) = 1 — S(z(t)) for 
tL < t < tff (Eriksson et al., 2012). 


5.3 Population dispersal 

We set up the initial population density at t = 50 kya as shown in the top 
left panel of Figure [9} The integration units are scaled following Tab. [l] such 
that A = 1.67 • 10 -3 ky _1 and c = 208 km 2 /ky (consistent with those used 
in (Young and Bettinger, 1995])). The solver is then run using time frame fC 


maps described above, down to 1 kya. The remaining plots (Figures [lofllj ) 
display the resulting population dispersal simulation on unsmoothed /C maps, 


regularized by (13). Using population parameters consistent with the liter¬ 


ature, the gross features of the late (50-60 kya) out-of-Africa dispersal of 


14 






















Homo sapiens are reproduced (Forster, 2004), e.g. the colonization of West¬ 
ern Europe by ~ 40 kya and that of South America before 14 kya. Using the 
solver with the same initial conditions but on a smooth NPP map, like the 
one shown in right panel of Figure [7J without the regularizer (13), yields the 
same wavefront propagation speed. 


6 Conclusions 

In this paper, we presented a novel semi-implicit Godunov scheme for the 
Fisher/KPP equation with a constant carrying capacity 1C, described in Sec- 
tion [3} In one dimension, the expected traveling wave ( [Fisher , 1937; Kol 


mogorov et ah, 1937) develops as shown in Figure [TJ In other tests, not 
shown here, we saw that almost any concentrated inital condition will de¬ 
velop similar waves: for example, two nearby peaks. Our scheme is on a 
rectangular grid, so in 2-D we need to ensure that in cylindrically symmetric 
situations we can control the errors due to the x — y asymmetry. The error 
plots, Figure [4j show that we can reduce the asymmetries somewhat by an 
alternating direction scheme. In any case, these errors are very small even 
for a relatively high CFL number. As in the one-dimensional situation, a 
traveling wave also develops as expected (Kolmogorov et ah, 1937): Figure[2] 
and Figure [3] show this development and compare the results to the 2-D sym¬ 
metric version of the Matlab function pdepe. Because the x discretization 
has truncation errors proportional to (Ax) 3 , after 0(|) time steps we should 
not be surprised to see the errors shown in Figure [3] turn up for very small h 
and behave roughly as 0(|(Ax) 3 ). 

In Secs. |4.1| we extended our procedure to handle an x-dependent 1C, 
specifically eq. ©• In order to regularize our solver against bad behavior 


when dealing with 1C maps inferred from real-world data, in Section 4.2 


we studied both regularization and the expected constant velocity of the 
traveling wave. Except for small transition regions entering or leaving a 
region of low carrying capacity, the velocity (y/2 in our scaling) is indeed 
constant. In Section |4~4[ we went further to develop a scheme for the situation 
with both x and t dependent carrying capacity 1C via eq.(16). 


Finally, this scheme has been applied to a prototypical case in population 
dynamics: the out-of-Africa dispersal of Homo sapiens. On the Mercator 
projected world map, by using vegetation net primal productivity as a proxy 
for carrying capacity, Figure [9] shows that by the regularization of 1C in space 
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via eq. (13) but interpolating in time yields stable and reasonable results. In 
fact, the results showing ancestor arrival in NE India at roughly 40 kya, then 
crossing the Bering Strait before 10 kya, and multiple routes into South Asia 
(Re yes-Centento et al.[ 2014) are very encouraging. Honesty requires that we 
admit our size (408km) 2 pixels do not resolve the two crossing points at Bab 
al Mendab and Sinai adequately. Additionally, a switch (see: Section 4.5) 
which would allow passage at the Bering Strait seems unnecessary due to the 
interesting coincidence that the hominin wave front reaches this passage at 
the beginning of the last ice age. If it were blocked previously, this would 
have had no effect. 

The core computations performed by our solver are independent tridiag¬ 
onal solutions, which can be easily parallelized to deal with larger grids. In 
order to further improve numerical performance, in the Appendix [XJ we dis¬ 
cussed a compressed storage scheme to integrate the Fisher/KPP equation 
on a projected world map (or any other irregularly-shaped domain). The 
alternating direction scheme discussed in Section |3.2| also works with this 
compressed storage. In the case of world dispersal discussed in Section [5j 
since about 71% of the earth’s surface is water, this compressed storage re¬ 
duces computational work by the same amount. 


Provenance 

For this paper, the simulations were run on either a Mac Mini, 2.4 GHz 
Intel Core Duo (Mac OS 10.6.8), or a MacBook laptop with the same proces¬ 
sor specifications but running OS 10.8.4. On the Mini, MatLab (7.10.0.499) 
R2010a was used, respectively (8.1.0.604) R2013a on the MacBook laptop. 
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A Map segmentation appendix 


In order to use our solver on a geographical map, it suffices to use a map 
outline, i.e., a rectangular grid with l’s in habitable regions, and 0’s in the 
water, as in Figure [T2j Since each Godunov direction step only involves a 
row, or a column, independently, we can set up the following indexing scheme. 
Each row i — 1... 50 in Figure 12 will have nysegs (i) of habitable segments, 
whose starting and ending positions are ystart_seg(k) and yend_seg(k), 
respectively, where k = 1... nysegs (i). Likewise, for j = 1... 100 columns 
each with nxsegs(j) also with start and end positions. Roughly 24% of the 
world map is land, i.e. habitable. As an example, Figure [12] shows row 26 has 
4 segments of varying size. Likewise, column 87 has 5 segments. A sample 
MatLab code in Appendix [A] illustrates the scheme for 1/2 step of x-direction 
updates, followed by a full step of y-direction updates, then again 1/2 step 
of x-direction segment updates. The alternating direction method, described 
in Section |T2[ just interchanges x ex y on alternate time steps. 

For the reader’s convenience, we include here a sample Matlab code of 
our Godunov-Strang-Yoshida scheme. 


% NY X direction updates for half-stepl 
locx = 0; 
for j=l:NY 

nsegs = nxsegs(j); 
for k=l:nsegs 

istart=xstart_seg(locx+k); iend=xend_seg(locx+k); 
ninseg=iend-istart+l; 

u0(l) = 0; u0(ninseg+2) =0; % boundary values 

u0(2:ninseg+l) = u(istart:iend,j); 

% eq. (15a) solution: 

ut = godunovstepl(ninseg+2,h,kcfl,u0,scl,sc2); 
u(istart:iend,j) = ut(2:ninseg+l); 

end 

locx = locx + nsegs; 

end 

% NX Y direction updates step2 
locy = 0; 
for i=l:NX 
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nsegs = nysegs(i); 
for k=l:nsegs 

jstart=ystart_seg(locy+k); jend=yend_seg(locy+k); 
ninseg=j end-jstart+1; 

ul(l) = 0; ul(ninseg+2) =0; °/ 0 boundary values 
ul(2:ninseg+l) = u(i,jstart:jend)’; 
cap(l) = 1; cap(ninseg+2) = 1; 
cap(2:ninseg+l) = kap(i,jstart:jend); 

% eq. (15c) solution: 

ut = godunovstep2(ninseg+2,h,kcfl,ul,scl,sc2,cap); 
u(i,jstart:jend) = ut(2:ninseg+l); 

end 

locy = locy + nsegs; 

end 

% Repeat godunovstepl, as above for eq. (15d) 

% locx = 0; 

l for j=l:NY 

% ETC 

% end 
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wave front at t=0 


wave front at t=30 




Figure 1: Left: Godunov vs. pdepe at the end of one time step t = h. Right: 
Godunov vs. pedpe at time t = 30. The time step h = 1/5 and Ax = 1/5. 


u(x,y) at t=5 u(x,y) at t=20 



Figure 2: Left: 2-D solution at t = 5 obtained with Godunov scheme 
Right: same, at t = 20. The time step was h = 1/10, Ax = 2/5 and 
N x = N y = 201. 
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Radial PDEPE vs. Godunov, t=20 



RMS error vs. h 



Figure 3: Left: 2-D Wave front at t = 20 compared to pdepe. Right: Rough 
RMS error as a function of step size, h, at t = 20. Errors e^MS and ^max ~ 
3.5 srms have approximately 0(h 2 ) behavior for h > 1/10, but increase if h 
is too small: overall estimate is 0(h 2 ) + 0(Ax 2 ) + 0(/Ax 3 ). 


error in u(t=20) 


error in u(t=20) 




Figure 4: Left: | u(x, y) — u(r) \ error prohle of the Godunov scheme compared 
to the pdepe solution at t = 20 symmetrizing x -H- y on alternate time steps: 
£rms = 7.3-10 -3 , £max = 2.9-10 -2 . Right: the same, but with fixed x and y: 
£rms = 7.4 • 10 -3 , 6max = 3.0 • 10 -2 . Step size is h = 1/10, N x = N y = 401. 
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log(1 + 10 u(x,y)) at t=37.5 
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Figure 5: Left: desert test density at t = 37.5: N x = 201, N y = 101, and 
h = 1/8. For x < —16 and x > 16, 1C = 1, while K, — f r — 0.01 in 
the |x| < 16 desert. The image shows u(x,y) shaded from 0 (black) to 1 
(white). The wave can be seen crossing into the upper /C = 1 region, while 
the population in the lower K = 1 region has already saturated. Right: half¬ 
height Xi/ 2 (t) of the wave front vs. time, compared with pdepe. Velocity V 
has a transient increase, then decrease, as the wave enters/exits the desert. 
Initial condition: a strip with ((x — Xi ) 2 ) 1 / 2 = 3 initial population starting 
at xi = —30. 


wave front x 1/2 vs. t, K=0.01 



regularizes {X=5) 



Figure 6: Regularizers for eq. (14). 
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Siberia raw (log), 50Kya 



Siberia smooth (log), 50Kya 



Figure 7: Left: Siberia rough 50kya NPP distribution. Note the K, = 0 hole 
at pixel north = 14, east = 38. Right: same sub-map but smoothed by the 
low pass filter (15). Scales are log(l + 10/C) with the original 0.3 < 1C < 
3 x 10 4 biomass units (gm/m 2 ) (Grieser et al., 
same colorbar. 


2006). Both plots use the 


0-1 sigmoids vs. straight line 



t 


Figure 8: Interpolation between /C maps at different times: sigmoids (18) for 


v = 1/2 and v = 1, straight line, and Eriksson’s /(/(/(£))) model (Eriksson 


et al. 2012). 
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Initial log(1 + 10 u(x,y)) att=50 Kya 



u -u at t= 1.00 Kya 

sm rg 1 



Figure 9: Color-coded log (1 + 10m) plot for out-of-Africa dispersal. Dark 
blue shows water , with JC = 0. Growth rate is A = L67 • 10 _3 ky _1 , and 
diffusion coefficient c = 208km 2 /ky (Young and Bettinger, 1995). Left: 
initial distribution at t = 50 kya. This color scale is used in all subsequent 
dispersal plots below. Right: differences between the regularization (13) of 
u vs. smoothing of /C(x, t) © at t = 1 kya. 


log(1 + 10 u(x,y)) at t=40 Kya 



log(1 + 10 u(x,y)) at t=30 Kya 



Figure 10: Out of Africa simulations by regularization (13). Left: population 
at t = 40 kya. Right: population at t = 30 kya. Color scale is from Fig. [9] 


(left). 
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log(1 + 10 u(x,y)) at t=11 Kya 


log(1 + 10 u(x,y)) at t=1.0 Kya 




Figure 11: Out of Africa simulations with regularization (13). Left: popu¬ 
lation at t = 11 kya. Right: population at t = 1 kya. Color scale is from 
Fig. @ (left). 


row 26, col 87 segmentation 



X 


Figure 12: Map with habitable regions set to 1 (light blue) and water regions 
set to 0 (dark blue) on a N x = 100, N y = 50 grid. The figure shows the 
segmentation for row 26 and column 87, with 4 and 5 segments, shown in 
cyan and yellow respectively. 
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